## Figure 3: Covid-19 Unemployment Effects Across Municipalities with Low and High Prejudice and Without and With Far-right Mayors

rm(list = ls())

setwd('/path/to/replication/')

library(data.table)
library(dplyr)
library(ggplot2)


if (!file.exists('./output/out_econ_het_vshare_con_trend_all.RData') &
    !file.exists('./output/out_econ_het_con_trend_all.RData') &
    !file.exists('./output/het_effects_ri.RData')) {
  source('./code/figure3reg.R')
}

# extract relevant estimates from regression models
out_econ_het_vshare_con_trend_all <- get(load('./output/out_econ_het_vshare_con_trend_all.RData'))
var <- (vcov(out_econ_het_vshare_con_trend_all)['covid2:losing_income_d','covid2:losing_income_d']) + (vcov(out_econ_het_vshare_con_trend_all)['covid2:losing_income_d:extr_right_d','covid2:losing_income_d:extr_right_d']) +
  (2*(vcov(out_econ_het_vshare_con_trend_all)['covid2:losing_income_d','covid2:losing_income_d:extr_right_d']))
sd <- sqrt(var)
estimate <- tidy(out_econ_het_vshare_con_trend_all) %>% filter(term==c('covid2:losing_income_d')) %>% select(estimate) +
  tidy(out_econ_het_vshare_con_trend_all) %>% filter(term==c('covid2:losing_income_d:extr_right_d')) %>% select(estimate)
conf.low <- estimate-1.96*sd
conf.high <- estimate+1.96*sd
estimates_econ_het_vshare_con_trend_all <- unname(cbind(estimate,conf.low, conf.high))
colnames(estimates_econ_het_vshare_con_trend_all) <- c('estimate', 'conf.low', 'conf.high')

out_econ_het_con_trend_all <- get(load('./output/out_econ_het_con_trend_all.RData'))
var <- (vcov(out_econ_het_con_trend_all)['covid2:losing_income_d','covid2:losing_income_d']) + (vcov(out_econ_het_con_trend_all)['covid2:losing_income_d:extr_right_mayor','covid2:losing_income_d:extr_right_mayor']) +
  (2*(vcov(out_econ_het_con_trend_all)['covid2:losing_income_d','covid2:losing_income_d:extr_right_mayor']))
sd <- sqrt(var)
estimate <- tidy(out_econ_het_con_trend_all) %>% filter(term==c('covid2:losing_income_d')) %>% select(estimate) +
  tidy(out_econ_het_con_trend_all) %>% filter(term==c('covid2:losing_income_d:extr_right_mayor')) %>% select(estimate)
conf.low <- estimate-1.96*sd
conf.high <- estimate+1.96*sd
estimates_econ_het_con_trend_all <- unname(cbind(estimate,conf.low, conf.high))
colnames(estimates_econ_het_con_trend_all) <- c('estimate', 'conf.low', 'conf.high')


# construct figure
out_econ_het <- rbind(tidy(out_econ_het_vshare_con_trend_all) %>% filter(term==c('covid2:losing_income_d')) %>% select(estimate, conf.low, conf.high),
                  estimates_econ_het_vshare_con_trend_all,
                  tidy(out_econ_het_con_trend_all) %>% filter(term==c('covid2:losing_income_d')) %>% select(estimate, conf.low, conf.high),
                  estimates_econ_het_con_trend_all)

types <- data.frame(group = c(paste('Low', 'prejudice', sep="\n"), paste('High', 'prejudice', sep="\n"),
                              paste('Without', 'far-right', 'mayor', sep="\n"), paste('With', 'far-right', 'mayor', sep="\n")),
                    color = c(rep(paste('Individual-level', 'mobilization', sep='\n'),2),
                              rep(paste('Local-level', 'mobilization', sep='\n'),2)))

out <- cbind(out_econ_het,types)
out$group <- factor(out$group, levels = c(paste('Low', 'prejudice', sep="\n"), paste('High', 'prejudice', sep="\n"),
                                          paste('Without', 'far-right', 'mayor', sep="\n"), paste('With', 'far-right', 'mayor', sep="\n")))
out$color <- factor(out$color, levels = c(paste('Individual-level', 'mobilization', sep='\n'),
                                          paste('Local-level', 'mobilization', sep='\n')))

out %>% ggplot() + aes(x = color, y = estimate, group = group, color=group) +
  geom_point(position=position_jitterdodge(seed=147)) +
  geom_hline(yintercept = 0, colour='#999999', linetype='longdash') + 
  geom_linerange(aes(ymin = conf.low, ymax = conf.high), position =position_jitterdodge(seed=147)) +
  scale_color_grey(start = 0, end = 0.8) +
  theme(
    legend.position="bottom", legend.title = element_blank(), legend.text = element_text(size = 12),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(), axis.line = element_line(colour = "black"
    ),axis.title.y = element_text(size = rel(1)), axis.text.x=element_text(size=12),axis.text.y=element_text(size=12)) +
  xlab(NULL)+
  ylab('Covid-19 Unemployment Effects on Asian Hate Crimes')
ggsave('./het_econ_effects.pdf', width=6, height=4.85)


## randomization inference tests on 'local-level mobilization'
# randomization inference test for effect in municipalities lead by far-right mayors
rr <- get(load('./output/het_effects_ri.RData'))
marginal_effect <- unname(unlist(rr$marginal_effect))
obs_marginal_effect <- unname(tidy(out_econ_het_con_trend_all) %>% filter(term==c('covid2:losing_income_d:extr_right_mayor')) %>% select(estimate))
p_value_marginal_effect <- (length(which(marginal_effect < -c(obs_marginal_effect[1,1])[[1]])) + length(which(marginal_effect > c(obs_marginal_effect[1,1])[[1]])))/length(marginal_effect)

# are the coefficients from Covid-19 unemployment effects across municipalities with/without far-right mayors distinguishable from each other?
beta1 <- tidy(out_econ_het_con_trend_all) %>% filter(term==c('covid2:losing_income_d')) %>% select(estimate)
beta2 <- estimate
var1 <- (tidy(out_econ_het_con_trend_all) %>% filter(term==c('covid2:losing_income_d')) %>% select(std.error))^2
var2 <- sd^2
Z <- (beta2-beta1)/(sqrt(var1+var2))

z <- unname(unlist(rr$z))
p_value_z <- (length(which(z > Z[1,1])))/length(z)

